Мы рассмотрим

  • Различные варианты анализа, применяющегося для тех случаев, когда зависимая перменная - счетная величина (целые неотрицательные числа)

Вы сможете

  • Объяснить особенности разных типов распределений, принадлежащих экспоненциальному семейству.
  • Построить пуасоновскую и квази-пуассоновскую линейную модель
  • Объяснить проблемы, связанные с избыточностью дисперсии в модели
  • Построить модель, основанную на отрицательном биномиальном распределении

Различные типы распределений

Распределение

То, что мы в быту привыкли называть распределением - это функция плотности вероятности.

Плотность вероятности - это функция, описывающая вероятность получения разных значений случайной величины

Нормальное распределение

\[f(y;\mu, \sigma)= \frac {1}{\sigma \sqrt{2 \pi}} e^{-\frac{(y-\mu)^2}{2\sigma^2}}\]

Два параметра (\(\mu\), \(\sigma\))

Среднее:   \(E(Y) = \mu\)

Дисперсия: \(var(Y) = \sigma^2\)

Пределы варьирования

\(-\infty \le Y \le +\infty\)

Распределение Пуассона

\[f(y;\mu)= \frac{\mu^y \times e{-\mu}}{y!}\]

Один параметр (\(\mu\))

Среднее:   \(E(Y) = \mu\)

Дисперсия: \(var(Y) = \mu\)

Важное свойство: При увеличении значения \(\mu\) увеличивается размах варьирования

Пределы варьирования

\(0 \le Y \le +\infty\),
\(Y\) целочисленные

Гамма-распределение

\(f(y; \mu, \nu) = \frac{1}{\Gamma(\nu)}\times (\frac{\nu}{\mu})^{\nu} \times y^{\nu-1} \times e^{\frac{y \times \nu}{\mu}}\)

Два параметра (\(\mu\), \(\nu\))

Среднее:   \(E(Y) = \mu\)

Дисперсия: \(var(Y) = \frac {\mu^2}{\nu}\)

Параметр \(\nu\) определяет степень избыточности дисперсии

Пределы варьирования

\(0 < Y \le +\infty\)
Внимание! \(Y\) строго больше 0

Отрицательное биномиальное распределение

\(f(y; k, \mu) = \frac{\Gamma(y + k)}{\Gamma(k) \times \Gamma(y+1)} \times (\frac{k}{\mu + k})^k \times (1 - \frac{k}{\mu + k})^y\)

Это смесь Пуассоновского и Гамма распределений: \(Y\) демонстрируют распределение Пуассона с \(\mu\), подчиняющимися Гамма-распределению.

Два параметра (\(\mu\), \(k\))

Среднее:   \(E(Y) = \mu\)
Дисперсия: \(var(Y) = \mu + \frac {\mu^2}{k}\)
Параметр \(k\) определяет степень избыточности дисперсии.
Важное свойство: Приближается к распр. Пуассона при очень больших \(k\).

Пределы варьирования

\(0 \le Y \le +\infty\),   \(Y\) целочисленные

Биномиальное распределение

\(f(y; N, \pi) = \frac{N!}{y! \times (N-y)!} \times \pi^y \times (1 - \pi)^{N-y}\)

Два параметра (\(N\), \(\pi\))

Среднее:    \(E(Y) = N \times \pi\)
Дисперсия: \(var(Y) = N \times \pi \times (1-\pi)\)
Параметр \(N\) определяет количество объектов в испытании
Парметр \(\pi\) - вероятность события (\(y = 1\))

Пределы варьирования

\(0 \le Y \le +\infty\),   \(Y\) целочисленные

Распределение зависимой переменной и линия регрессии

Связующая функция (Link function)

Наиболее распространенные (канонические) связующие функции

Характер величин Распределение Связующая функция (link function)
Непрерывные величины, потенциально варьирующие в пределах \(-\infty , + \infty\) Гауссовское (Нормальное) identity \(X\beta = \mu\)
Бинарные величины (1; 0), или количество (доля) объектов, относящихся к одному из двух классов Биномиальное распределение logit \(X\beta = ln(\frac{\mu}{1 - \mu})\)
Счетные величины (0, 1, 2, 3…) Распределение Пуассона или Отрицательное биномиальное распределение log \(X\beta = ln(\mu)\)

NB! Есть и другие связующие функции

Связующая функция (Link function)

Связующая функция (Link function)

Связующая функция (Link function)

Модели, основанные на распределении Пуассона и отрицательном биномиальном распределении

Способствуют ли взрослые мидии притоку молоди?

Даные взяты из работы Khaitov, 2013

juv_ad <- read.table("data/mussel_juv_ad.csv", sep=";", header =T)
head(juv_ad, 12)

##    Year Bank Sample Juv Adult
## 1  1997 vor4      1  90    42
## 2  1997 vor4      2 134    25
## 3  1997 vor4      3 166    27
## 4  1997 vor4      4 168    48
## 5  1997 vor4      5 102    25
## 6  1997 vor4      6  69    57
## 7  1998 vor4      1 410    21
## 8  1998 vor4      2  73    76
## 9  1998 vor4      3 347     0
## 10 1998 vor4      5 402     5
## 11 1998 vor4      6 158    28
## 12 1999 vor4      1 282     0

В этих данных ожидается проблема автокорреляции остатков

Вопрос:

Как можно учесть в модели многолетний характер сбора материала?

Построим простую линейную модель

M1 <- glm(Juv ~ Adult * factor(Year), data = juv_ad)
drop1(M1, test = "F")
## Single term deletions
## 
## Model:
## Juv ~ Adult * factor(Year)
##                    Df Deviance  AIC F value Pr(>F)  
## <none>                  256799 1024                 
## Adult:factor(Year) 14   361458 1026    1.72  0.076 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Можно убрать взаимодействие

M2 <- lm(Juv ~ Adult  +  factor(Year), data = juv_ad)
library(car)
Anova(M2)
## Anova Table (Type II tests)
## 
## Response: Juv
##               Sum Sq Df F value  Pr(>F)    
## Adult          63886  1    12.9 0.00059 ***
## factor(Year) 1357772 14    19.6 < 2e-16 ***
## Residuals     361458 73                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Посмотрим на предсказания этой модели

  • Модель предсказывает, что взрослые негативно влияют на обилие молоди.
  • Модель предсказывает отрицательные значения!

Диагностика модели

M2_diag <- fortify(M2)

ggplot(M2_diag, aes(x = .fitted, y = .stdresid)) + geom_point() +
  geom_hline(yintercept = 0) + geom_smooth(se = FALSE)

  • Явные признаки гетероскедастичности!

Два способа решения проблем с моделью

  1. Провести логарифмирование зависимой переменной и построить модель для логарифмированных величин.
  2. Построить модель, основанную на распределении Пуассона.

Модель, основанная на распределении Пуассона

M3 <- glm(Juv ~ Adult * factor(Year), 
          data = juv_ad, family = "poisson")
library(car)
Anova(M3)
## Analysis of Deviance Table (Type II tests)
## 
## Response: Juv
##                    LR Chisq Df Pr(>Chisq)    
## Adult                   272  1     <2e-16 ***
## factor(Year)           5031 14     <2e-16 ***
## Adult:factor(Year)      439 14     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Диагностика модели

M3_diag <- data.frame(.fitted = predict(M3),
                      .pears_resid = residuals(M3, type = "pearson"))

Рассеяние остатков выглядит лучше!

Избыточность дисперсии (Overdispersion)

В Пуассоновской регрессии мы моделируем изменение распределения Пуассона в зависимости от каких-то предикторов.

В распределении Пуассона \(E(Y) = \mu\) и \(var(Y) = \mu\)

Если в имеющихся данных \(var(Y) > \mu\), то нарушается условие применимости пуассоновской регрессии.

Избыточность дисперсии (Overdispersion)

Первый способ оценки избыточности дисперсии

Resid_M3 <- resid(M3, type = "pearson") # Пирсоновские остатки

N <- nrow(juv_ad) # Объем выборки

p <- length(coef(M3))   # Число параметров в модели

df <- (N - p) # число степенейсвободы

fi <- sum(Resid_M3^2) /df  #Величина fi показывает во сколько раз в среднем sigma > mu для данной модели

fi
## [1] 16.8

Дисперсия в 16.756 раза больше среднего!

Избыточность дисперсии (Overdispersion)

Второй способ оценки избыточности дисперсии

library(qcc)
qcc.overdispersion.test(juv_ad$Juv, type = "poisson")
##                    
## Overdispersion test Obs.Var/Theor.Var Statistic p-value
##        poisson data               103      9079       0

Избыточность дисперсии (Overdispersion)

Очень маленькие стандартные ошибки (и все очень достоверно) - это явный признак избыточности дисперсии

summary(M3)
## 
## Call:
## glm(formula = Juv ~ Adult * factor(Year), family = "poisson", 
##     data = juv_ad)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -10.779   -2.370   -0.074    2.127    7.858  
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             5.14952    0.11551   44.58  < 2e-16 ***
## Adult                  -0.00955    0.00305   -3.13  0.00173 ** 
## factor(Year)1998        0.89375    0.12042    7.42  1.2e-13 ***
## factor(Year)1999        0.46774    0.12279    3.81  0.00014 ***
## factor(Year)2000       -0.67585    0.20462   -3.30  0.00096 ***
## factor(Year)2001        0.09944    0.23147    0.43  0.66748    
## factor(Year)2002       -0.46104    0.25068   -1.84  0.06589 .  
## factor(Year)2003        0.99365    0.12334    8.06  7.9e-16 ***
## factor(Year)2004        0.34742    0.13293    2.61  0.00896 ** 
## factor(Year)2005        0.51396    0.13149    3.91  9.3e-05 ***
## factor(Year)2006        0.40032    0.12997    3.08  0.00207 ** 
## factor(Year)2007       -1.51444    0.21138   -7.16  7.8e-13 ***
## factor(Year)2008        0.89713    0.12124    7.40  1.4e-13 ***
## factor(Year)2009        0.47749    0.12401    3.85  0.00012 ***
## factor(Year)2010        0.31445    0.13054    2.41  0.01600 *  
## factor(Year)2011        1.05365    0.12258    8.60  < 2e-16 ***
## Adult:factor(Year)1998 -0.01150    0.00336   -3.42  0.00063 ***
## Adult:factor(Year)1999  0.01496    0.01048    1.43  0.15355    
## Adult:factor(Year)2000  0.02677    0.00580    4.62  3.9e-06 ***
## Adult:factor(Year)2001 -0.00481    0.00490   -0.98  0.32698    
## Adult:factor(Year)2002  0.00534    0.00533    1.00  0.31643    
## Adult:factor(Year)2003 -0.02341    0.00411   -5.69  1.2e-08 ***
## Adult:factor(Year)2004 -0.01631    0.00591   -2.76  0.00575 ** 
## Adult:factor(Year)2005 -0.00882    0.00381   -2.31  0.02062 *  
## Adult:factor(Year)2006 -0.04550    0.00760   -5.99  2.1e-09 ***
## Adult:factor(Year)2007  0.05039    0.00944    5.34  9.4e-08 ***
## Adult:factor(Year)2008  0.00569    0.00349    1.63  0.10312    
## Adult:factor(Year)2009 -0.00642    0.00508   -1.26  0.20637    
## Adult:factor(Year)2010  0.00445    0.00398    1.12  0.26394    
## Adult:factor(Year)2011  0.01949    0.00351    5.55  2.8e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 7992.3  on 88  degrees of freedom
## Residual deviance: 1008.2  on 59  degrees of freedom
## AIC: 1695
## 
## Number of Fisher Scoring iterations: 4

Источники избыточности дисперсии

  1. Мнимая избыточность дисперсии
    • Наличие отскакивающих значений
    • Как следствие пропущенных ковариат или взаимодействий предикторов
    • Наличие внутригрупповых корреляций (нарушение независимости выборок)
    • Нелинейный характер взаимосвязи между ковариатами и зависимой переменной
    • Неверно подобранная связывающая функция
    • Количество нулей больше, чем предсказывает распределение Пуассона (Zero inflation)
  2. Истинная избыточность дисперсии, как следствие природы данных.

Как бороться с избыточностью дисперсии

Если избыточнсть дисперсии мнимая, то ее надо устранить, введя в модель соответствующие поправки.

Если избыточность дисперсии истинная, то необходима более серьезная коррекция модели.

Два пути решения проблемы при истинной избыточности дисперсии

  1. Построить квази-пуассоновскую модель
  2. Построить модель, основанную на отрицательном биномиальном распределении

Квази-пуассоновские модели

Отличие от пуассоновсой модели заключается лишь в том, что в квази-пуассоновских моделях вводится поправка для связи дисперсии и матожидания.

В этой модели матожидание \(E(Y) = \mu\) и дисперсия \(var(Y) = \phi \times \mu\)


Величина \(\phi\) показывает во сколько раз дисперсия превышает матожидание.

\[\phi = \frac{var(Y)}{\mu}=\frac {\frac{\sum{(\epsilon_i)^2}}{N - p}} {\mu} = \frac{\sum{(\epsilon_{pearson})^2}}{N - p}\]

  • Модель, по сути, остается той же, что и пуассоновская, но изменяются стандартные ошибки оценок параметров, они домножаются на \(\sqrt{\phi}\)
  • Для квази-пуассоновских моделй не определена функция максимального правдоподобия и, следовательно, нельзя вычислить AIC

Квази-пуассоновская модель

M4 <- glm(Juv ~ Adult * factor(Year), data = juv_ad, family = "quasipoisson")
Anova(M4)
## Analysis of Deviance Table (Type II tests)
## 
## Response: Juv
##                    LR Chisq Df Pr(>Chisq)    
## Adult                  16.2  1   0.000056 ***
## factor(Year)          300.3 14    < 2e-16 ***
## Adult:factor(Year)     26.2 14      0.024 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Важно: Распределение разности девианс лишь приблизительно описывается распределением \(\chi^2\). Поэтому уровни значимости близкие к 0.05 нельзя рассматривать как надежные.

Уровень значимости для взаимодействия в данной модели близок к 0.05. Можно подумать об упрощении модели!

Квази-пуассоновская модель

Упрощенная модель

M4a <- glm(Juv ~ Adult  + factor(Year), data = juv_ad, family = "quasipoisson")
Anova(M4a)
## Analysis of Deviance Table (Type II tests)
## 
## Response: Juv
##              LR Chisq Df Pr(>Chisq)    
## Adult            13.8  1     0.0002 ***
## factor(Year)    254.9 14     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Предсказания упрощенной модели

  • Биологический вывод: взрослые мидии препятствуют пополнению молодью

Модель, основанная на отрицательном биномиальном распределении

library(MASS)
M5 <- glm.nb(Juv ~ Adult*factor(Year) , data = juv_ad, link = "log")
Anova(M5)
## Analysis of Deviance Table (Type II tests)
## 
## Response: Juv
##                    LR Chisq Df Pr(>Chisq)    
## Adult                    27  1 0.00000016 ***
## factor(Year)            370 14    < 2e-16 ***
## Adult:factor(Year)       34 14     0.0023 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Уровень значимости для взаимодействия заметно ниже 0.05. Взаимодействие факторов отбросить нельзя!

Задание

Проверьте на избыточность дисперсии модель, основанную на отрицательном биномиальном распеделении

Решение

Resid_M5 <- resid(M5, type = "pearson") # Пирсоновские остатки
N <- nrow(juv_ad) # Объем выборки
p <- length(coef(M5)) +1   # Число параметров в модели
df <- (N - p) # число степенейсвободы
fi <- sum(Resid_M5^2) /df  #Величина fi показывает
# во сколько раз в среднем sigma > mu для данной модели
fi
## [1] 1.47

Диагностика модели

M5_diag <- data.frame(.fitted = predict(M5), 
                      .pears_resid = residuals(M5, type = "pearson"))

ggplot(M5_diag, aes(x=.fitted, y = .pears_resid)) + geom_point() +
  geom_hline(yintercept = 0) + geom_smooth(se = F)

Задание

Визуализируйте предсказания модели, основанной на отрицательном биномиальном распределении

Визуализируем предсказание модели

  • Биологический вывод: В разные годы харакер взаимосвязи молоди и взрослых мидий может быть разным

Код для графика

MyData <- unique(juv_ad)


MyData$Predicted <- predict(M5, newdata = MyData, type = "response")
ggplot(MyData, aes(x = Adult, y = Predicted, group = Year)) + 
    geom_line(color = "blue") + geom_hline(yintercept = 0) + 
    ylab("Ожидаемое количество молоди") + geom_point(data = juv_ad, 
    aes(x = Adult, y = Juv, color = factor(Year))) + 
    facet_wrap(~Year, ncol = 5) + guides(color = FALSE)

Таким образом

  1. Модели, основанные на неподходящем типе распределения, могут давать нереалистичные предсказанные значения.
  2. В зависимости от того, как сконструирована модель, можно получить результаты, которые позволят сформулировать принципиально разные биологические выводы.

Выбор оптимальной модели

Какие факторы определяют супружескую неверность?

data(Affairs, package = "AER")
af <- Affairs

Оригинальная работа: Fair, R.C. (1978). A Theory of Extramarital Affairs. Journal of Political Economy, 86, 45–61.

affairs - Количество внебрачных свзяей за последний год
gender - пол
age - возраст
yearsmarried - сколько ле в браке
children - наличие детей
religiousness - уровеь религиозности
education - уровень образования
rating - самооценка ощущений от брака

Задание:

  1. Постройте оптимальную модель, описывающую зависимость количества внебрачных связей в зависимости от пола, времени, проведенного в браке, наличия детей, уровня религиозности и уровня образованности.
  2. Проверьте валидность данной модели

Решение

Mod <- glm(affairs ~ gender * yearsmarried * children * religiousness + 
    education, data = af, family = "poisson")

summary(Mod)
## 
## Call:
## glm(formula = affairs ~ gender * yearsmarried * children * religiousness + 
##     education, family = "poisson", data = af)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
##  -3.27   -1.78   -1.26   -0.60    6.74  
## 
## Coefficients:
##                                                   Estimate Std. Error
## (Intercept)                                         0.6798     0.5380
## gendermale                                         -1.2262     0.7071
## yearsmarried                                        0.2747     0.0666
## childrenyes                                         0.7479     0.6325
## religiousness                                      -0.5268     0.1915
## education                                          -0.0274     0.0149
## gendermale:yearsmarried                             0.0325     0.0796
## gendermale:childrenyes                              1.7851     0.8714
## yearsmarried:childrenyes                           -0.2238     0.0740
## gendermale:religiousness                            0.5167     0.2605
## yearsmarried:religiousness                         -0.0183     0.0241
## childrenyes:religiousness                           0.1421     0.2405
## gendermale:yearsmarried:childrenyes                -0.0833     0.0902
## gendermale:yearsmarried:religiousness              -0.0282     0.0288
## gendermale:childrenyes:religiousness               -0.7454     0.3216
## yearsmarried:childrenyes:religiousness              0.0216     0.0266
## gendermale:yearsmarried:childrenyes:religiousness   0.0499     0.0324
##                                                   z value Pr(>|z|)
## (Intercept)                                          1.26   0.2064
## gendermale                                          -1.73   0.0829
## yearsmarried                                         4.13 0.000037
## childrenyes                                          1.18   0.2371
## religiousness                                       -2.75   0.0059
## education                                           -1.84   0.0662
## gendermale:yearsmarried                              0.41   0.6829
## gendermale:childrenyes                               2.05   0.0405
## yearsmarried:childrenyes                            -3.03   0.0025
## gendermale:religiousness                             1.98   0.0473
## yearsmarried:religiousness                          -0.76   0.4473
## childrenyes:religiousness                            0.59   0.5544
## gendermale:yearsmarried:childrenyes                 -0.92   0.3555
## gendermale:yearsmarried:religiousness               -0.98   0.3274
## gendermale:childrenyes:religiousness                -2.32   0.0205
## yearsmarried:childrenyes:religiousness               0.81   0.4174
## gendermale:yearsmarried:childrenyes:religiousness    1.54   0.1242
##                                                      
## (Intercept)                                          
## gendermale                                        .  
## yearsmarried                                      ***
## childrenyes                                          
## religiousness                                     ** 
## education                                         .  
## gendermale:yearsmarried                              
## gendermale:childrenyes                            *  
## yearsmarried:childrenyes                          ** 
## gendermale:religiousness                          *  
## yearsmarried:religiousness                           
## childrenyes:religiousness                            
## gendermale:yearsmarried:childrenyes                  
## gendermale:yearsmarried:religiousness                
## gendermale:childrenyes:religiousness              *  
## yearsmarried:childrenyes:religiousness               
## gendermale:yearsmarried:childrenyes:religiousness    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2925.5  on 600  degrees of freedom
## Residual deviance: 2518.8  on 584  degrees of freedom
## AIC: 3047
## 
## Number of Fisher Scoring iterations: 7

Проверка на избыточность дисперсии

# Проверка на Overdispersion
Resid_Mod <- resid(Mod, type = "pearson") 
N <- nrow(af) 
p <- length(coef(Mod)) 
df <- (N - p) 
fi <- sum(Resid_Mod^2) /df  
fi
## [1] 6.44

Строим квази-пуассоновскую модель

Mod1 <- glm(affairs ~ gender * yearsmarried * children * religiousness * 
    education, data = af, family = "quasipoisson")

summary(Mod1)
## 
## Call:
## glm(formula = affairs ~ gender * yearsmarried * children * religiousness * 
##     education, family = "quasipoisson", data = af)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -3.356  -1.699  -1.197  -0.518   6.833  
## 
## Coefficients:
##                                                              Estimate
## (Intercept)                                                  -0.72294
## gendermale                                                    6.32757
## yearsmarried                                                 -0.13878
## childrenyes                                                  21.26596
## religiousness                                                 2.68502
## education                                                     0.10930
## gendermale:yearsmarried                                       1.87507
## gendermale:childrenyes                                      -27.44754
## yearsmarried:childrenyes                                     -1.24552
## gendermale:religiousness                                     -3.35695
## yearsmarried:religiousness                                   -0.01597
## childrenyes:religiousness                                   -10.10729
## gendermale:education                                         -0.52222
## yearsmarried:education                                        0.01782
## childrenyes:education                                        -1.44076
## religiousness:education                                      -0.23255
## gendermale:yearsmarried:childrenyes                          -0.86718
## gendermale:yearsmarried:religiousness                        -0.47529
## gendermale:childrenyes:religiousness                         11.26910
## yearsmarried:childrenyes:religiousness                        0.56754
## gendermale:yearsmarried:education                            -0.09418
## gendermale:childrenyes:education                              1.97936
## yearsmarried:childrenyes:education                            0.08050
## gendermale:religiousness:education                            0.27615
## yearsmarried:religiousness:education                          0.00383
## childrenyes:religiousness:education                           0.71028
## gendermale:yearsmarried:childrenyes:religiousness            -0.00469
## gendermale:yearsmarried:childrenyes:education                 0.01680
## gendermale:yearsmarried:religiousness:education               0.01994
## gendermale:childrenyes:religiousness:education               -0.82013
## yearsmarried:childrenyes:religiousness:education             -0.04129
## gendermale:yearsmarried:childrenyes:religiousness:education   0.01511
##                                                             Std. Error
## (Intercept)                                                   11.55629
## gendermale                                                    15.47884
## yearsmarried                                                   1.28616
## childrenyes                                                   14.96782
## religiousness                                                  4.61169
## education                                                      0.77941
## gendermale:yearsmarried                                        2.14783
## gendermale:childrenyes                                        19.12367
## yearsmarried:childrenyes                                       1.45496
## gendermale:religiousness                                       5.84741
## yearsmarried:religiousness                                     0.47583
## childrenyes:religiousness                                      5.74652
## gendermale:education                                           1.00639
## yearsmarried:education                                         0.08122
## childrenyes:education                                          1.01421
## religiousness:education                                        0.31845
## gendermale:yearsmarried:childrenyes                            2.31968
## gendermale:yearsmarried:religiousness                          0.70212
## gendermale:childrenyes:religiousness                           7.13962
## yearsmarried:childrenyes:religiousness                         0.53301
## gendermale:yearsmarried:education                              0.12702
## gendermale:childrenyes:education                               1.24743
## yearsmarried:childrenyes:education                             0.09344
## gendermale:religiousness:education                             0.38721
## yearsmarried:religiousness:education                           0.03019
## childrenyes:religiousness:education                            0.39169
## gendermale:yearsmarried:childrenyes:religiousness              0.76524
## gendermale:yearsmarried:childrenyes:education                  0.13881
## gendermale:yearsmarried:religiousness:education                0.04196
## gendermale:childrenyes:religiousness:education                 0.46874
## yearsmarried:childrenyes:religiousness:education               0.03419
## gendermale:yearsmarried:childrenyes:religiousness:education    0.04622
##                                                             t value
## (Intercept)                                                   -0.06
## gendermale                                                     0.41
## yearsmarried                                                  -0.11
## childrenyes                                                    1.42
## religiousness                                                  0.58
## education                                                      0.14
## gendermale:yearsmarried                                        0.87
## gendermale:childrenyes                                        -1.44
## yearsmarried:childrenyes                                      -0.86
## gendermale:religiousness                                      -0.57
## yearsmarried:religiousness                                    -0.03
## childrenyes:religiousness                                     -1.76
## gendermale:education                                          -0.52
## yearsmarried:education                                         0.22
## childrenyes:education                                         -1.42
## religiousness:education                                       -0.73
## gendermale:yearsmarried:childrenyes                           -0.37
## gendermale:yearsmarried:religiousness                         -0.68
## gendermale:childrenyes:religiousness                           1.58
## yearsmarried:childrenyes:religiousness                         1.06
## gendermale:yearsmarried:education                             -0.74
## gendermale:childrenyes:education                               1.59
## yearsmarried:childrenyes:education                             0.86
## gendermale:religiousness:education                             0.71
## yearsmarried:religiousness:education                           0.13
## childrenyes:religiousness:education                            1.81
## gendermale:yearsmarried:childrenyes:religiousness             -0.01
## gendermale:yearsmarried:childrenyes:education                  0.12
## gendermale:yearsmarried:religiousness:education                0.48
## gendermale:childrenyes:religiousness:education                -1.75
## yearsmarried:childrenyes:religiousness:education              -1.21
## gendermale:yearsmarried:childrenyes:religiousness:education    0.33
##                                                             Pr(>|t|)
## (Intercept)                                                    0.950
## gendermale                                                     0.683
## yearsmarried                                                   0.914
## childrenyes                                                    0.156
## religiousness                                                  0.561
## education                                                      0.889
## gendermale:yearsmarried                                        0.383
## gendermale:childrenyes                                         0.152
## yearsmarried:childrenyes                                       0.392
## gendermale:religiousness                                       0.566
## yearsmarried:religiousness                                     0.973
## childrenyes:religiousness                                      0.079
## gendermale:education                                           0.604
## yearsmarried:education                                         0.826
## childrenyes:education                                          0.156
## religiousness:education                                        0.466
## gendermale:yearsmarried:childrenyes                            0.709
## gendermale:yearsmarried:religiousness                          0.499
## gendermale:childrenyes:religiousness                           0.115
## yearsmarried:childrenyes:religiousness                         0.287
## gendermale:yearsmarried:education                              0.459
## gendermale:childrenyes:education                               0.113
## yearsmarried:childrenyes:education                             0.389
## gendermale:religiousness:education                             0.476
## yearsmarried:religiousness:education                           0.899
## childrenyes:religiousness:education                            0.070
## gendermale:yearsmarried:childrenyes:religiousness              0.995
## gendermale:yearsmarried:childrenyes:education                  0.904
## gendermale:yearsmarried:religiousness:education                0.635
## gendermale:childrenyes:religiousness:education                 0.081
## yearsmarried:childrenyes:religiousness:education               0.228
## gendermale:yearsmarried:childrenyes:religiousness:education    0.744
##                                                              
## (Intercept)                                                  
## gendermale                                                   
## yearsmarried                                                 
## childrenyes                                                  
## religiousness                                                
## education                                                    
## gendermale:yearsmarried                                      
## gendermale:childrenyes                                       
## yearsmarried:childrenyes                                     
## gendermale:religiousness                                     
## yearsmarried:religiousness                                   
## childrenyes:religiousness                                   .
## gendermale:education                                         
## yearsmarried:education                                       
## childrenyes:education                                        
## religiousness:education                                      
## gendermale:yearsmarried:childrenyes                          
## gendermale:yearsmarried:religiousness                        
## gendermale:childrenyes:religiousness                         
## yearsmarried:childrenyes:religiousness                       
## gendermale:yearsmarried:education                            
## gendermale:childrenyes:education                             
## yearsmarried:childrenyes:education                           
## gendermale:religiousness:education                           
## yearsmarried:religiousness:education                         
## childrenyes:religiousness:education                         .
## gendermale:yearsmarried:childrenyes:religiousness            
## gendermale:yearsmarried:childrenyes:education                
## gendermale:yearsmarried:religiousness:education              
## gendermale:childrenyes:religiousness:education              .
## yearsmarried:childrenyes:religiousness:education             
## gendermale:yearsmarried:childrenyes:religiousness:education  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 6.52)
## 
##     Null deviance: 2925.5  on 600  degrees of freedom
## Residual deviance: 2378.1  on 569  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 9

Подбираем оптимальную модель

step(Mod1) #Для квази-пуассоновских моделей эта функция работать
# не будет, так как не определен AIC

Строим модель, основанную на отрицательном биномиальном распределении

Mod_nb <- glm.nb(affairs ~ gender * yearsmarried * children * 
    religiousness + education, data = af)
Anova(Mod_nb)
## Analysis of Deviance Table (Type II tests)
## 
## Response: affairs
##                                            LR Chisq Df Pr(>Chisq)    
## gender                                         0.12  1     0.7282    
## yearsmarried                                  17.87  1   0.000024 ***
## children                                       1.24  1     0.2647    
## religiousness                                 15.06  1     0.0001 ***
## education                                      1.93  1     0.1645    
## gender:yearsmarried                            0.22  1     0.6366    
## gender:children                                0.02  1     0.8975    
## yearsmarried:children                          7.98  1     0.0047 ** 
## gender:religiousness                           0.01  1     0.9179    
## yearsmarried:religiousness                     0.28  1     0.5965    
## children:religiousness                         0.02  1     0.8758    
## gender:yearsmarried:children                   0.82  1     0.3662    
## gender:yearsmarried:religiousness              0.80  1     0.3715    
## gender:children:religiousness                  2.25  1     0.1337    
## yearsmarried:children:religiousness            1.78  1     0.1825    
## gender:yearsmarried:children:religiousness     0.06  1     0.8032    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Подбираем оптимальную модель

step(Mod_nb)

Сморим на результаты оптимальной модели

Mod_nb_final <- glm.nb(formula = affairs ~ yearsmarried + 
    children + religiousness + yearsmarried:children, 
    data = af, init.theta = 0.1346363532, link = log)

summary(Mod_nb_final)
## 
## Call:
## glm.nb(formula = affairs ~ yearsmarried + children + religiousness + 
##     yearsmarried:children, data = af, init.theta = 0.1346359871, 
##     link = log)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.090  -0.824  -0.700  -0.372   1.910  
## 
## Coefficients:
##                          Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)               -0.4047     0.4309   -0.94    0.3476    
## yearsmarried               0.2978     0.0631    4.72 0.0000024 ***
## childrenyes                1.3847     0.4621    3.00    0.0027 ** 
## religiousness             -0.4171     0.1064   -3.92 0.0000878 ***
## yearsmarried:childrenyes  -0.2232     0.0690   -3.24    0.0012 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.135) family taken to be 1)
## 
##     Null deviance: 375.76  on 600  degrees of freedom
## Residual deviance: 335.57  on 596  degrees of freedom
## AIC: 1477
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.1346 
##           Std. Err.:  0.0148 
## 
##  2 x log-likelihood:  -1465.2630

Проводим диагностику оптимальной модели

Mod_nb_test <- data.frame(.fitted = predict(Mod_nb_final), 
                          .pears_resid = residuals(Mod_nb, type = "pearson"))

ggplot(Mod_nb_test, aes(x=.fitted, y = .pears_resid)) + geom_point() +
  geom_smooth(se = FALSE)

Проверим на избыточность дисперсии

Resid_Mod <- resid(Mod_nb_final, type = "pearson") 
N <- nrow(af) 
p <- length(coef(Mod_nb_final)) 
df <- (N - p) 
fi <- sum(Resid_Mod^2) /df  
fi
## [1] 0.723

Визуализируем предсказание модели

Код для графика

MyData <- expand.grid(yearsmarried = seq(min(af$yearsmarried), 
    max(af$yearsmarried)), children = c("yes", "no"), 
    religiousness = seq(min(af$religiousness), max(af$religiousness)))
MyData$Predicted <- predict(Mod_nb_final, newdata = MyData, 
    type = "response")

ggplot(MyData, aes(x = yearsmarried, y = Predicted, 
    color = religiousness)) + geom_line(aes(group = religiousness), 
    size = 1.5) + facet_grid(~children, labeller = label_both) + 
    scale_color_gradient(low = "yellow", high = "red") + 
    geom_point(data = af, aes(x = yearsmarried, y = affairs), 
        position = position_jitter(width = 1, height = 1))

Summary

  • В случае счетных зависимых перменных (неотрицательных целочисленных величин) применяются модели, основанные на распределении Пуассона или отрицаетльном биномиальном распределении.
  • Важным ограничивающим условием применения этих моделей является отсутствие избыточности дисперсии.
  • Избыточность дисперсии может быть истинной и мнимой.
  • При истинной избыточности дисперсии модель можно скорректировать, поcтроив квази-пуассоновскую модель (вводятся поправки для ошибок оценок коэффициентов модели).
  • Другой подход - построение моделей, основанных на отрицательном биномиальном распределении.

Что почитать

  • Кабаков Р.И. R в действии. Анализ и визуализация данных на языке R. М.: ДМК Пресс, 2014.
  • Zuur, A.F. et al. 2009. Mixed effects models and extensions in ecology with R. - Statistics for biology and health. Springer, New York, NY.